% This function computes term premia. For 
% TP_version == 1 we use the defition of term premium in Rudebusch and
% Swanson. For TP_version = 2 we use the definition in Dai and Singleton
% (2002) where the convexity term is included in the term premium

function model = getTermPremia_loadings(model,setupFilter)


% Solve for expected short rates
var0   = model.by0(1,1);
varx   = model.byx(1,:);
varxx  = model.byxx(1,:,:);
varxxx = model.byxxx(1,:,:,:);
Max_maturity = max(setupFilter.maturities);
[r0,rx,rxx,rxxx]= projectionStepExpectedControlClosedForm(var0,varx,varxx,varxxx,...
    max(setupFilter.maturities),model,setupFilter.setupProj);
nx           = size(model.hx,1);


%% The Term Premium
% Allocating memory for the output
TP0   = zeros(Max_maturity,1);
TPx   = zeros(Max_maturity,nx);
TPxx  = zeros(Max_maturity,nx^2);
TPxxx = zeros(Max_maturity,nx^3);
% Term premia = R_long - E_t*average(R_t+i-1,i=1:n))
% For the nominal term premia
cum_sum0           = zeros(Max_maturity,1);
cum_sumx           = zeros(Max_maturity,nx);
cum_sumxx          = zeros(Max_maturity,nx^2);
cum_sumxxx         = zeros(Max_maturity,nx^3);
% The first term is special:
cum_sum0(1,1)      = var0;
cum_sumx(1,:)      = varx;
cum_sumxx(1,:)     = reshape(varxx,1,nx^2);
cum_sumxxx(1,:)    = reshape(varxxx,1,nx^3);
for i=2:Max_maturity
    cum_sum0(i,1)       = cum_sum0(i-1,1)    + r0(i-1,1);
    cum_sumx(i,:)       = cum_sumx(i-1,:)    + rx(i-1,:);
    cum_sumxx(i,:)      = cum_sumxx(i-1,:)   + reshape(rxx(i-1,:,:),1,nx^2);
    cum_sumxxx(i,:)     = cum_sumxxx(i-1,:)  + reshape(rxxx(i-1,:,:,:),1,nx^3);
end
for i=1:Max_maturity
    TP0(i,:)       = model.by0(i,:)       - 1/i*cum_sum0(i,:);
    TPx(i,:)       = model.byx(i,:)       - 1/i*cum_sumx(i,:);
    TPxx(i,:)      = model.byxx(i,:)      - 1/i*cum_sumxx(i,:);
    TPxxx(i,:)     = model.byxxx(i,:)     - 1/i*cum_sumxxx(i,:);
end

%% The output
model.TP0    = TP0;
model.TPx    = TPx;
model.TPxx   = TPxx;
model.TPxxx  = TPxxx;

end